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We present a semiclassical theory of linear and nonlinear optical response of graphene. The em¬ 
phasis is placed on the nonlinear optical response of graphene from the standpoint of the underlying 
chiral symmetry. The Bloch quasiparticles in low energy limit, around the degeneracy points are 
dominantly chiral. It is shown for the first time that this chiral behavior in conjunction with scale 
invariance in graphene around the Dirac points results in the strong nonlinear optical response. 

Explicit expressions for the linear and nonlinear conductivity tensors are derived based on Semicon¬ 
ductor Bloch Equations (SBEs). The linear terms agree with the result of Kubo formulation. The 
three main additive mechanisms contribute in the nonlinear optical response of graphene: pure in¬ 
traband, pure interband and the interplay between them. For each contribution, an explicit response 
function is derived. The Kerr-type nonlinearity of graphene is then studied and it is demonstrated 
that its Kerr nonlinear coefficient is several orders of magnitude higher than that of many other 
known semiconductors. In addition, the nonlinear refractive index of graphene can also be tuned 
and enhanced by applying a gate voltage. 


I. INTRODUCTION 

Graphene is a two dimensional arrangement of the car¬ 
bon atoms sitting in a hexagonal lattice, a seemingly sim¬ 
ple lattice structure that nonetheless underlies the special 
transport and optical properties [l|. The band struc¬ 
ture of graphene differs substantially from other con¬ 
densed matter systems. The effective Hamiltonian de¬ 
scribes pseudo-relativistic quasiparticles obeying (2 -|- I)- 
dimensional Dirac equation. In the context of QED, the 
electronic excitations introduced by such dynamics can 
be considered as massless chiral fermions 0. 

Graphene exhibits a variety of peculiar properties 
that are manifestations of the special symmetries of its 
crystalline structure and relativistic energy spectrum of 
charged carriers. Symmetries entail several unconven¬ 
tional properties such as the existence of a topologi¬ 
cally protected zero-energy state. Berry phase, anoma¬ 
lous quantum Hall effect and Zitterbewegung (‘trembling 
motion’) [H, H, 0|- It is counterintuitive that intrinsic 
graphene has a finite conductivity, in the order of the 
Hall conductivity at zero temperature and zero car¬ 
rier concentration. The current operator does not com¬ 
mute with the Hamiltonian and therefore, graphene can¬ 
not sustain the current. This intrisic disorder leads to 
a finite conductivity All these odd properties can be 
linked to the chiral behavior of the carriers. In graphene 
the pseudospin is locked parallel or antiparallel to the di¬ 
rection along which the electron propagates and so the 
quasiparticles possess the property of chirality 

The optical response of graphene is also expected to 
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be influenced by the chiral nature of the carriers and the 
scale invariance of the band structure in low energy limit. 
However, despite its importance, a theoretical study on 
the unconventional optical response of graphene is still 
lacking. The optical response of graphene in the linear 
regime has been investigated theoretically and experi¬ 
mentally iSIl- Graphene as a scale invariant two di¬ 
mensional chiral electronic system exhibits universal op¬ 
tical response [l^ . A simple analysis based on linear 
response theory shows that an isolated graphene sheet 
can absorb about 2.3% of the normally incident opti¬ 
cal field, which is indeed a huge number for a mono- 
layer atomic structure. The nonlinear optical response 
of graphene has been a topic of intensive research in the 
recent years Treatment of the nonlinear opti¬ 

cal response of graphene in the framework of quasiclassi- 
cal Boltzman equation predicts strong nonlinearity in the 
terahertz range of frequenc y, n eglecting pair productions 
and interband transitions [l3l| . This part of nonlinear¬ 
ity is mainly due to the geometrical properties of band 
structure rather than its topological aspects [l^. The 
calculation of optical response of graphene in time do¬ 
main has been carried out in Ref. |l5l |. Wrighte et al. 
have performed Fourier analysis of the dirac equation to 
obtain the optical response of the system for a given inci¬ 
dent field [14| . All existing time domain methods proceed 
primarily at the level of the wavefunction, rather than at 
the level of the density matrix and thereby suffer from 
the computational cost and the difficulty resulted from 
the inclusion of relaxation processes due to impurities 
and emission. The nonlinear optical resposen in gapped 
graphene, where the low-energy single-particle spectrum 
is modeled by massive Dirac equation is discussed in [2lj . 

The treatment of optical response of graphene, in a 
simplistic manner, proceeds from single particle approx¬ 
imation and the equation of motion for density matrix. 
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The formal approach to calculate the optical response 
of semiconductors, excited by a electromagnetic field, is 
based on perturbative expansion of the density matrix 
and taking all possible transitions into account [l^, [2^ . 
Even at such a reasonably simplistic treatment, the eval¬ 
uation of nonlinear response coefficients is numerically 
difficult and does not provide any intuitive insight. To 
circumvent these difficulties. Semiconductor Bloch Equa¬ 
tions (SBEs) for graphene are employed. 

SBEs [15 have extensively been used to calculate the 
nonlinear optical response of semiconductors [1^ . [27j| . In 
Ref. [15 the derivation of SBEs for graphene beyond 
the Dirac cone approximation has been discussed. In 
Ref. [15 the general treatment of the SBEs for graphene 
including electron-electron interactions and exciton ef¬ 
fects is presented. Using SBEs, the problem of interac¬ 
tion can be treated in a semiclassical manner leading to 
numerically amenable expressions for arbitrary orders of 
interaction. SBEs introduce an effective dipole in the re¬ 
ciprocal space revealing peculiarities of graphene in terms 
of its optical response. We have shown that, in the ab¬ 
sence of spin interactions, the aforementioned dipole is 
singular at the high symmetry points of the reciprocal 
space leading to a strong nonlinear optical response. 

In this paper we study the optical response of graphene 
based on a semicassical theory. We will adopt an ap¬ 
proach that treats the electrons dynamics in the pres¬ 
ence of a moderate intensity electromagnetic field based 
on SBEs. We demonstrate that the higher order non¬ 
linear response of graphene, possesses a singularity due 
to the topological properties of the band structure and 
the chiral nature of the charged carriers. To remove this 
singularity the inclusion of the spin-orbit interaction in 
a phenomenological level is proposed. 

The paper is organized as follows. In Sec. [Ill we present 
the Hamiltonian of graphene within the tight binding ap¬ 
proximation. In Sec. mil we address the question of how 
the chirality of carriers affects the optical response of 
graphene and its dependence on the fundamental group 
symmetries of the problem. In Sec. m the equations of 
motion for the single particle density matrix are formu¬ 
lated and the SBEs are derived in the context of Dirac 
equation. In Sec. El we propose an iterative approach 
to solve the effective optical Bloch equations. Secs. EH 
and IVIII present the derivation of linear and third order 
nonlinear conductivities. The numerical results and dis¬ 
cussions on the importance of graphene as a strong non¬ 
linear material are given in Sec. IVIIIl We will show that 
nonlinear effects in graphene are substantially stronger 
than those of other known semiconductors. In Sec. in 
we summarize our results. 


II. GRAPHENE HAMILTONIAN AND 
EQUATIONS OF MOTION 

Graphene has a honeycomb crystal lattice with two 
lattice points per elementary cell. They belong to two 


sublattices A and B where the nearest neighbours of the 
sites of one of them are sites belonging to the other sub¬ 
lattice (a bipartite lattice). In Fig. |l(a)] atoms in A and B 
sublattices are shown by blue and red balls respectively. 
The Bravais lattice is triangular with the lattice vectors 
given by Q, 



As shown in Fig. |l(b)| the reciprocal lattice is also hexag¬ 
onal with rhomboidal unit cell formed by two vectors 

The high symmetry crystallographic points are pre¬ 
sented in Fig. |l(b)| Throughout this paper the graphene 
monolayer (laying on the xy-plane) interacts with a plane 
wave illuminating graphene in the perpendicular direc¬ 
tion. This assumption allows us to use the electric dipole 
approximation in which the effect of magnetic field is 
excluded. Actually this approximation is quite accurate 
for an ideal graphene sheet wherein electrons are strongly 
bounded and their off-plane dynamics is negligible. The 
electric field can have an arbitrary time variation con¬ 
taining different harmonics. The dynamical properties of 
the positively charged ions that constitute the host lat¬ 
tice of the crystal will be neglected in our formulations. 
The system Hamiltonian for a single graphene sheet in¬ 
teracting with a classical electromagnetic field within the 
single particle approximation is 

H = Ho + Hi (3) 

Where Hq governs the dynamics of the electrons with 
mass M in the presence of the periodic lattice potential 
U(r) 


Ho = J d^r^^ir) + U(r)| ^(r) (4) 

The interaction Hamiltonian in long-wavelength limit (or 
normal incidence) is rigorously obtained in the velocity- 
gauge by replacing p by p -I- eA in Hp where A is the 
associated vector magnetic potential [^. In the case of 
graphene it seems very simple to use this electrodynam¬ 
ics substitution, however it can be shown that neither 
the calculations are efficient, nor it reveals some interest¬ 
ing physical properties. The interaction problem can be 
recast into the length gauge 


Hi = eEj„(t) 


d^r4'^(r)r4'(r) 


( 5 ) 


A common problem with the perturbation theory for 
solids in the length gauge is the difficult treatment of the 

H ition operator r in view of the extended Bloch states 
. This is further detailed in the next section. It will 
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FIG. 1. The graphene lattice (a) and its reciprocal lattice (b) , The Dirac points K' S'lid K — 3 ^^ j are shown 

in the figure .The schematic plot of dispersion relation around the Dirac point (c) 


be shown that both Hq and Hj can be treated in the 
framework of the Tight-Binding (TB) regime. Following 
the TB approximation the field operator is expanded in 
terms of 2 pz orbital wavefunction (j)(y) 

^(r) = X! X] - R.A)ak 

k Ra 

+ -X ^ ^ exp(ik.RB)(/i(r - RB)Sk 
k Rb 

= + ^B(r) ( 6 ) 

The summation run over all sublattices coordinates de¬ 
noted by 'R.a/b- The operators hk and 6 k are fermionic 
annihilation operators on A and B sublattices respec¬ 
tively. The TB hopping parameters can be found by 
fitting the results of the first-principles electronic struc¬ 
ture calculations with the experimental results. The sim¬ 
plest TB Hamiltonian contains hopping to the nearest- 
neighbor sites. The resulted Hamiltonian is easy to di¬ 
agonalized in the second quantized operators hk and 6k- 
Inserting (IHl) into (|1]) and taking into account the nearest- 
neighbor hopping K ~ —2.97eV [s^ along 82 and ^3 
bonds (shown in FiglT]) the TB Hamiltonian in the mo¬ 
mentum space is 


Ho = '^Eo -I- 6j(.6k) 

k 

+ kX + /*(k)6j^ak 

k 


( 7 ) 


Eq is the result of the hopping processes within the sub¬ 
lattices. The first term on the right hand side of the 
Eq. d?!) is symmetrically diagonal and does not affect 
the quasi-particles dynamics. The function /(k) carries 
the symmetry properties of the graphene lattice and is 
given by 
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f(^)='^exp{ik.S^) ( 8 ) 

i=l 


The TB Hamiltonian Hq becomes diagonal in the con¬ 
duction and valence basis 

ac = ^ + e+^-^/^Sk) (9) 

a. = ^ (e-“>^/ 2 ak - ( 10 ) 

where /(k) = |/(k)| exp(tak). This yields the corre¬ 
sponding TB-based band structure 

E^^/^ = Eo±n\f{k)\ (11) 

It is easy to show that the two bands cross at K and K' 
and the Hamiltonian displays the interesting physics of 
the Dirac fermions in the proximity of the conical points. 
These special crossing points are topologically protected 
and no band gap can be opened in the presence of any 
kind of perturbation preserving time reversal and inver¬ 
sion symmetries. For the case of intrinsic graphene sheet 
the Fermi surface coincides with the energy at the coni¬ 
cal points. For a reasonably doped graphene sheet, just 
the dynamics of the quasiparticles around the conical 
points play the major role in the electronic properties 
of graphene. In this scenario, graphene can be viewed as 
a vanishing-gap semiconductor. Some other crystalline 
structures such as HgTe are also known to be gapless 
semiconductors (3lj |. but what makes graphene unique is 
also the helical nature of the quasiparticles In the 
next sections it will be shown that this helicity plays a 
decisive role in nonlinear optical response of graphene 
which is the focus of this paper. In the vicinity of the 
conical points /(k) can be linearly expanded in terms of 
k components as 

/(K -H k) « (12) 


where (/9k is the angle of vector k with respect to the kx 
axis shown in Fig. |l(c) Linearization of the Hamiltonian 
around the conical points yields a massless Dirac quasi¬ 
particle whose dispersion relation is E\^ = Ehvpk. It is 
noted that vp = —3aK/2h is the Fermi velocity which 
is around c/300. Within the band structure picture the 
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effective Hamiltonian is 

^ ~ hvpk ~ (13) 

k 

where ^kc and ^kii are the conduction and valence anni¬ 
hilation operators in the upper and lower energy bands, 
respectively. As long as the inter-valley scattering is im¬ 
probable the local behavior of the Hamiltonian around 
K and K' is independent. For mathematical convenience 
the matrix representations are used in the derivation of 
the equations of motion. In fact, the system can be ad¬ 
equately described in atomistic language. Assume that 
two-component spinors (l O) and (O l) are adopted 
for A and B states respectively. The resulting time de- 
pendeted Dirac equation describing low energy excitation 
around one of the conical points is written as 

Hk = hvp^ ■ CT 

= ■Hkl'k (14) 

ot 

where 4^k is a two-component spinor in the momen¬ 
tum space and a is made from the Pauli matrices a = 
xOg. yUy + zUz arose form the two sublattices. The 
optical response of graphene can be studied in a physi¬ 
cally transparent way using the equation of motion for 
the density matrix. In Sec. llVI we examine the nonlinear 
optical response of graphene based on the time evolution 
of the density matrix. Liouville’s equation governs the 
time progress of the density matrix and it yields coupled 
differential equations. 

d ^ - 

*fi^/Ok = [H, Pk] = [Hu, Pk] + [Hi, pk] (15) 

The first term on the right hand side of Eq. (II3 is the 
regular dynamical phase variation. In the next section 
it will be shown that the second term is closely related 
to the Berry connection and topological properties of the 
band structure. This is the point where chiral nature of 
the carriers and strong optical response of graphene tie 
up. In the next section, an intuitive clue will be provided 
in a general formalism of two-band systems. 


III. TWO-BAND SYSTEMS AND THE ROLE OE 
CHIRALITY 

To illustrate the impact of the chiral nature of the 
charged carriers on the optical response of graphene and 
to explore the uniqueness of the graphene in terms of its 
strong nonlinear interaction with light, the mathematical 
description of chirality for a general two-level systems is 
presented. We also address applicability of the reduced 
TB basis to describe the matrix elements of the interac¬ 
tion Hamiltonian in the length gauge and its connection 
with chirality of the charged carriers. In the last part of 
this section, it is shown that, our arguments are general 


enough and they are independent of the approximation 
existing in TB calculations. Chirality and its influence 
on the optical response root in the discrete symmetries 
existing in the crystalline structure of graphene and in¬ 
clusion of the many body effects and the other higher 
order interaction terms do not alter the general conclu¬ 
sion. 


A. Two-band systems 


A prototypical two-band system might be described by 
the Hamiltonian expanded in terms of the Pauli matrices 
as 

H = eo(k)/ -I- £(k)M(k).(T (16) 

where eo(k) and e(k) > 0 are real functions of Bloch 
wavenumber k. The three-dimensional vector operator a 
is made from the Pauli matrices and / is a 2 x 2 iden¬ 
tity matrix. The vector 'u(k) is a three dimensional unit 
vector and can be represented in terms of the spherical 
angle variables a and /3 

U:j(k) = sin/3(k) cosa(k) (17) 

'Uy(k) = sin/?(k) sina(k) (18) 

'u^(k) = cos/3(k) (19) 


By successive applications of the rotation operator 
= exp(— ifi.CT^) (h and (j) are the axis and the 
ang le of the rotation respectively) along the Euler axis 
[3^ the eigenvectors correspond to two energy eigenval¬ 
ues E±(k) = eo(k) ± e(k) of the Hamiltonian can be 
obtained 


|k,t) 

|k,;) 


cos f I j e 
sin j e +“/2 

sin (f) 

— cos g +*“/2 


( 20 ) 

( 21 ) 


Eor a d-dimensional electronic system the Bloch mo¬ 
mentum k can be represented by its magnitude k 
and d — 1 angle variables in the spherical coordinates, 
■ ■ ■ , 7 d-i}- For the sake of brevity all angle vari¬ 
ables are conveniently called 7 . For the particular case 
of graphene the system is two-dimensional and only the 
azimuthal angle (^k is needed to determine the direction 
of the Bloch momentum in the reciprocal space. Suppose 
a prototype Hamiltonian in which u is a function of angle 
variables only, i.e. u = u^y), with no dependency on k, 
this is known as a general ehiral system [l2|. For such a 
system a and /3 appearing in Eqs. (120(1 and ((211) only de¬ 
pend on the angle variables. Equivalently the pseudospin 
is determined by the direction of the momentum. 

According to Eq. (fT4ll it is obvious that the low energy 
Hamiltonian in graphene describes a scale invariance chi¬ 
ral electronic system. The chiral symmetry of the carri¬ 
ers is not restricted to the TB model but stems from 






5 


the honeycomb translational symmetry of the crystalline 
structure. 


B. Position operator: The role of chirality 


To examine the importance of chirality in the opti¬ 
cal response of graphene, we now turn our attention 
to the calculation of matrix elements of the interaction 
Hamiltonian in the length gauge. As mentioned ear¬ 
lier, the calculation of matrix elements of the position 
operator in different Bloch states is challenging and it 
has caused some controversies [13, . For the general 

case of extended Bloch states with spatial dependency 
of 'l'„k(r) = (r|n,k) = exp('ik.r)M„k(r) where u„k is the 
periodic part of the wavefunction, the matrix elements 
of the pos ition operator are related to Berry connection 
tensor ]33j | . It is also shown [s^ that [f, pk] appeared on 
the left hand side of Eq. (USD can be expressed as 

[i*: Pk] — [Ak, Pk] 

where Ak,nm = —*(u„,k|Vk|um,k) is the Berry connec¬ 
tion tensor. Actually, the most difficulty in calculation 
of the matrix elements of the position operator is due to 
difficulty in the calculations of Berry connection. In order 
to compute the Berry connection tensor rigorously, the 
full machinery of Density Function Th eory and Wannier 
interpolation scheme should be used [S^. Introducing 
maximally localized Wannier basis functions provides a 
numerically feasible scheme to evaluate the matrix ele¬ 
ments. It is straightforward to show that in the proper 
gauge that the basis functions are expanded around the 
atomic centers Tq, as 

^n.k(i’) = ^ C'„,„k exp [ik.(R -h r^)] (/ia(r - R - r^) 

R,,q 

the right hand side of Eq. dSlD reads 


[f, Pk] = -iVkPk + [Ck: Pk] (23) 


where Ck.nm = -* Y.a C'a,nkVkC'a,mk E closely related to 
the Berry connection tensor. Localization of basis func¬ 
tions is the basic principle that underlines this approxi¬ 
mation. Fortunately, for the case of graphene the basis 
functions are fairly well localized and this approxima¬ 
tion works well. It is worth mentioning that both terms 
appearing on the right hand side of Eq. (l23ll are gauge 
dependent, but the overall expression is independent of 
gauge and the specific choice of basis functions. For the 
particular case of two band chiral systems described in 
the previous section, the Berry connection exhibits singu¬ 
lar behavior at the degeneracy points. Energy eigenstates 
only depend on the angular variables q^’s. Therefore the 
gradient operator acting on the angular functions will be 


VkC, 


Q;,kn 


1 

k 



i 


^ 


a 


OiMn 


(24) 


where khi{'^) and ji are respectively, the Riemann scale 
function and the unit vector associated with the angle 
variable ■ji. The most interesting property of the chiral 
systems in their optical response originate from the 1 /k 
dependence. The appearance of I/A: term in the Liouville 
equation is the main difference between graphene and an 
ordinary semiconductor material. It is shown in the next 
section that this term acts like a dipole in the reciprocal 
space, playing a significant role in the graphene’s nonlin¬ 
ear optical response. 


C. Chirality and symmetry 

Basically, low energy excitations that capture the uni¬ 
versal characteristics of the system are highly influenced 
by symmetries. Let us focus on the effective Hamiltonian 
governing electrons dynamics around the Fermi energy 
level including band renormalizations due to electron- 
electron interactions. In this paper we will not plunge 
into the Landau theory and just symmetry considerations 
are discussed here. As long as the two Dirac points can 
be treated as independent entities, first order expansion 
of the Hamiltonian around the conical points reads 

~ ^ + m{k)az (25) 

where i and j run over x and y. The coefficients A^-’s are 
the elements of a 2 x 2 constant real matrix. The mass 
term m(k) can be expanded as m(k) = mo -I- k^mx + 
kyuriy. Based on the mathematical description of the chi¬ 
rality elucidated above, it is easy to show that the nec¬ 
essary and sufficient condition to have a chiral system in 
low energy limit is mo = 0, which implies gapless state. 
It can be shown that the Dirac fermions (gapless prop¬ 
erty) are topologically protected Two main sym¬ 

metries characterize the hexagonal lattice with identical 
atoms on A and B sites. The first one is time reversal 
which exists in the absence of the magnetic interactions 
and complex hopping. This symmetry is present regard¬ 
less of space symmetry properties of atomic potential. 
The second symmetry, inversion, is induced by the mirror 
symmetry of the atomic arrangement. Clearly this sym¬ 
metry is present in the isotropic hopping scenario. It is 
easy to show that time reversal and inversion symmetries 
separately establish a relationship between the Hamilto¬ 
nians for different valleys and only if both symmetries 
are preserved ttiq must necessarily vanish and the Dirac 
nodes would be locally stable. Any kind of perturbation 
satisfying these symmetries cannot open a gap as long as 
the Dirac points do not meet up. Moreover, two inequiv¬ 
alent Dirac points carry vortices with opposite signs and 
the nodes with opposite vortices cannot be removed by 
themselves. For large enough perturbation, however, two 
different Dirac points may meet each other and annihi¬ 
late the vortices and open up a gap. It can be shown 
that in the presence of C^h symmetry the gapless Dirac 
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points are globally stable [s^. Taking all these symme¬ 
tries into account the low energy Hamiltonian describes 
chiral quasiparticles around the Dirac points. 

In the subsequent sections it will be shown that this 
chiral behavior leads to a singularity in higher order op¬ 
tical response of graphene. To remove this singular be¬ 
havior while keeping the symmetries intact a suitable 
band renormalizations can be used. In 2005, Kane and 
Mele showed that at sufficiently low energy an isolated 
graphene exhibits a quantum spin Hall effect with an en¬ 
ergy gap induced by spin-orbit interaction [35| . Although 
this band gap is very small, graphene as a critical elec¬ 
tronic state can be strongly affected by this perturbation. 
Taking all contributions into account, the spin-orbit cou¬ 
pling Hamitonian is 

Hso — ^so^z'^z^z T {^x'^z^y dyTo^a,) (^6) 

where di , fi and Si are Pauli matrices acting on pseu¬ 
dospin, valley index and electron spin respectively. The 
coefficients Ago and Xr are microscopic spin-orbit cou¬ 
pling constant and Rashba coefficient, (as a result of 
breaking mirror symmetry) respectively. The spin-orbit 
coupling factor Ago can be affected by curvature of the 
graphene sheet. The reported value for this coefficient for 
the ideal case of flat defect-free graphene is Ago ~ l/reV 
[ 3 ^. For all practical structures Xr is much smaller than 
Ago and the resulting energy gap is 2{Ago — Xr). This 
small gap can remove the singularity in the optical re¬ 
sponse of graphene which will be discussed in the next 
sections. 


IV. SEMICONDUCTOR BLOCH EQUATIONS 
FOR GRAPHENE 

As mentioned in the preceding sections, within the sin¬ 
gle particle approximation, the density matrix obeys dy¬ 
namical equations in Schrddinger’s picture. The applied 
field drives the distribution out of equilibrium leading to 
nonvanishing induced current. The dynamical equations 
on the density matrix can be recast into the form of the 
Semiconductor Bloch Equations (SBEs). Based on SBEs 
the dynamics is governed by a quasiclassical theory with 
quantum fluctuations superimposed. The quantum cor¬ 
rections to the classical dynamics will be converted to the 
well known problem of light-atom interaction . 


A. Equations of motion 

We proceed from Eq. (1^^ which offers a gauge inde¬ 
pendent relation and thus we are at liberty to choose any 
kind of gauge making the mathematical structure sim¬ 
pler. Working in the sublattice (pseudospin) basis and 
making use of Eq. (USD and Eq. (1^ gives 

[d, pk] -I- feE. VkPk (27) 


Due to the smallness of the band gap induced by spin- 
orbit coupling, the dispersion properties of the charged 
carriers would barely deviate from massless relativistic 
dynamics and can be safely neglected around the Dirac 
points. However, we have to acknowledge the impor¬ 
tance of this effect in the nonlinear optical response of 
graphene. The phenomenological inclusion of spin-orbit 
coupling as well as scattering due to imperfections will 
be discussed at the end of this section. The 2x2 pseu¬ 
dospin density matrix pk can be expanded in terms of 
Pauli matrices 


/3k = nwl + mk- a 


(28) 


On substituting Eq. (1^ into Eq. (EZD, one obtains de¬ 
coupled equations for charge density Uk and pseudospin 
density @ 


dt 

i9mk 

~dr 


-E.VkUk 

n 

g 

2 vf (k X TOk) -I- —E.Vkffi-k 


(29) 

(30) 


The right hand side of Eq. (1501) is analogous to spin pro¬ 
cession in a magnetic field. The same can be set for the 
pseudospin in the psudomagnetic field acting in the re¬ 
ciprocal space 0. This equation encodes a wealth of 
information about the optical response of graphene in¬ 
cluding linear and nonlinear response in noninteracting 
regime. Owing to the linear dispersion relation around 
the Dirac points, current operator has only paramagnetic 
component 


ff ea-Hk c. 

and the current density becomes 
J = (^k) = Tr {xx + yy)-^ mk (32) 

k 

Having derived the equations of motion in the sublattice 
basis, now, we can switch to the energy diagonal basis. 
To avoid confusion, we use to denote the matrix 
representation of the operators in the valence and con¬ 
duction basis. In the energy diagonal basis (l O) and 

(O l)^, stand for the upper and the lower energy lev¬ 
els respectively. In the energy diagonal basis the density 
matrix and the current operator become: 

Pk =Ink +rfik. (kaz + (fkay - zax'^ (33) 

Jk=- evF (kcr^ -I- ipkCTy'j (34) 

Where k and (pk are shown in Fig. |l(c)| In the thermal 
equilibrium, before switching on the incident field, the 
density distribution obeys fermion statistics 

{iUkc)o = fim) , {iUkx)o = fi-£(k)) (35) 
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Where subscript 0 denotes equilibrium state and £’(k) = 
hvfk is the upper energy level. The distribution f{E) is 
the fermionic distribution function 


l + exp(|^) 

where ^ and T are, respectively, the chemical potential 
associated with the fermi energy level Ef and the tem¬ 
perature. 


B. Dynamics of population difference and 
polarization 

In the presence of electromagnetic held, the current 
operator acquires a hnite expectation value. The parti- 

I 


aAr(k, t) 

dt 


fE.VkAr(k,t) 


dV{k,t) 

dt 


tE.VkiP(k,t) 


where <I)(k, t) is the effective dipole causing interband 
transitions in the reciprocal space and fl(k) is the Rabi 
frequency associated with interband transitions 

d>(k,t) = |^ (39) 

fl(k) = = 2vFk (40) 

h 

The coupled equations given in (l38ll are called semi¬ 
conductor Bloch equations (SBEs) for graphene. These 
equations must be solved simultaneously, subject to the 
initial condition imposed by the fermion distribution be¬ 
fore turning on the held. 

Mik-oo) = fi£{k))-f{-£{k)) 
iP(k,-oo) = 0 

The left side of the SBEs are essentially similar to the 
semiclassical Boltzman’s transport equation. This part 
of dynamics is responsible for intraband transitions for 
a pure graphene, neglecting the effect of collisions and 
imperfections. A simple way of incorporating the effect 
of collision into the theory is to use a complex frequency 
in the spectral domain. The right side of SBEs appear to 
resemble the problem of two level atomic transition in the 
presence of an applied electric held. However, the dipole 
that causes transition has been replaced by $(k, t) in the 
reciprocal space. The chiral nature of the charged carri¬ 
ers leaves its hngerprint on the appearance of 1/fc in the 
effective dipole expression. This singular behavior roots 
in 1/fc dependence of Ck in th® energy diagonal basis. As 


cle current, can be divided into two distinct parts. The 
hrst part is the current resulting from disturbing the dis¬ 
tribution of the charged carriers residing on the upper 
and lower energy levels and the second contribution is 
due to interference between them. The former is intra¬ 
band and the latter is the interband current. Following 
this statement, it will be shown that the optical response 
of a general two level system depends on the population 
difference A/'(k) and polarization Vik) 

■^(k, t) = (CLCkc) - = 2k.m (36) 

^(k, t) = = -z.m + (37) 

Taking M and V as dynamical varibales and using 
Eq. (1301) . we obtain the equations of motion for the pop¬ 
ulation difference and the polarization. 


= —2<l)(k, t)Im{7^(k, t)} 

in{k)r{k,t) -h |$(k,t)Ar(k,t) 


(38) 


discussed earlier, this singularity can be resolved by the 
phenomenological inclusion of spin-orbit coupling. Due 
to smallness of this effect, the induced mass can modify 
the effective dipole expression as 


$(k,t) 


e E ■ (/'k 
h v/fc2 + (5fc)2 


(41) 


where hvpSk = (A^o — Xu). The theory developed here 
serves as the starting point to analyze the optical re¬ 
sponse of graphene for an arbitrary order of interaction. 


V. SOLUTION TO THE SEMICONDUCTOR 
BLOCH EQUATIONS 

As discussed, the SBEs in their original form de¬ 
scribe the quasiclassical transport and interband excita¬ 
tion problems simultaneously. To convert the dynamical 
equations into a more convenient form, we proceed to de¬ 
couple the transport and interband evolutions. It is noted 
that, neglecting the right side of SBE’s, the Boltzman- 
type transport equation introduces a moving frame in the 
reciprocal space responsible for intraband evolution. It is 
worth noting that for a moderate applied field strength, 
the time evolution due to this moving frame can be con¬ 
sidered as an adiabatic evolution in comparison to the in¬ 
terband one. Following this adiabatic argument, we now 
proceed to decouple the theory by introducing a moving 
frame in the reciprocal space. Assume that ko (t) defines 












the primed frame as 


Mww-coordinate system. 


e /■* 

ko{t) = --J Edt (42) 

then k = ko + k'. The equation of motion in the primed 
frame is merely that of a simple two level problem 


= -2$(k'+ ko,t)Im{P(k'+ ko,f)} 
= +in(k^ + kQ)7^(k^ + kQ, t) 

+f$(k' + ko,t)AA(k' + ko,t) 

(43) 

These equations resemble th e op tical Bloch equations for 
a generic two level problem [24| . In compliance with the 
standard notation of Bloch equations, we introduce u>k, 
Uk, Uk as 


cIA/’(k'+ ko, t) 

m 

^ drik' + ko,t) 
dt 


Wk>(t) = Af {k'+ko(t),t) (44) 

Uk'(t) = 2Re{P(k'+ ko(t),t)} (45) 

'Ck'(i) =-2Im{P(k'+ ko(t),t)} (46) 

Wk is population inversion in the moving frame. In 
Eqs. (H5]| and (HSl) . the factor of two and the minus sign 
are used to comply with convention, ^k and ojk are also 
defined for mathematical convenience: 

^k'(t)=<k(k' + ko(t),t) (47) 

uJk'(t) = fl(k' + ko(t)) (48) 


Wk = -7i (wk - Wk'^) + ^k-Ck (49) 

Uk = Wk-Ck - 72“k (50) 

Wk = -WkWk - 72Wk - CkWk (51) 

Where ‘dot’ denotes time derivative. The function is 
the population difference at equilibrium: 

= f(£(k)) - f(-£(k)) (52) 

For a weak pump field, the inversion Wk tends to relax 
to w^'^. The coherent terms, on the other hand, are the 
oscilatory functions of the field. To proceed further, the 
current response in the reciprocal space must be identi¬ 
fied. According to Eq. (1551) together with Eq. (IMl) the 
induced current is 


J = —2evF + <^k<,2k)- 


ruk 


= evF ^ l^-Wk-kok + Wk-kfliPk (53) 

k 


Therefore the equation of motion describing the time evo¬ 
lution of Uk provides enough information to model the in¬ 
terband response of graphene. Neglecting the time vari¬ 
ations of LOk in adiabatic approximation: 

Wk+272Uk+(wk + 72 ) Wk = -l2^kWk-ikWk-^kWk (54) 

Since is much larger than 7 I , we can drop 7 |wk to 
obtain the result: 


The functions ^k and Wk are the analytical functions of 
the exciting field. According to Eq. gm, ^k is the equiv¬ 
alent dipole in the moving frame. This dipole explicitly 
depends on the exciting field in the numerator of Eq. dm) 
and higher order nonlinear terms also exist due to motion 
of the primed frame. In adiabatic approximation the dy¬ 
namical equations (not the excitations) are not directly 
affected by the moving frame and therefore the time vari¬ 
ations of Wk can be neglected as long as the pump wave 
intensity is not so large that multiphoton excitations take 
place. 

The equations (l4^ provide an adequate description 
of resonant interband optical processes under conditions 
where relaxation processes can be neglected. Like the 
other two-level problems, the consideration of the relax¬ 
ation processes, typically due to spontaneous emission, 
is of central importance [ 2 ^. We assume that the upper 
energy level decays to the lower energy level at a rate of 
7 i. We also assume that the coherence terms decay with 
the dephasing rate of 72 . Due to the strong k dependence 
of the transition dipole in the reciprocal space, the decay 
rates are expected to be k dependent. However, in the 
simplest approach, these parameters are assumed to be 
constant. The coupled Bloch equations can be converted 
to the well known optical Bloch equations in two level 
approximation [^ . From now on, we drop the prime in 


Wk + 272Wk + WkWk = - 72 ^kWk - CkWk - CkWk (55) 

Eq. (1551) describes a driven damped harmonic oscilla¬ 
tor problem. The master equation (1551) in conjunction 
with Eq. (l4^ describe all linear and nonlinear proper¬ 
ties of graphene. The origin of the nonlinear interband 
response in the moving frame lies in the fact that the 
coupling to the optical field depends parametrically on 
the inversion Wk- Inversion is driven by the field stength 
^k as described by Eq. (l4^ which leads to a pure inter¬ 
band nonlinearity. This set of equations can be solved 
iteratively. Expanding Vk and Wk into the powers of the 
exciting field, i.e. ^ki gives a infinite series that contains 
the odd powers for Vk and even powers of the field for 
Wk- 

00 

n;k = <" + E (56) 

n—1 

00 

= E (57) 

n—1 

In the rest of the paper, the n’th order expansion terms 
w^^ and are defined via 







9 


In addition to the pure interband multiphoton pro¬ 
cess described above, a part of nonlinearity originates 
from the quasiclassical transport or intraband transi¬ 
tions. As will be clarified further in section IVIIl the 
frequency mixing effects in graphene arise from the pure 
intraband, pure interband and interband-intraband tran¬ 
sitions. The pure intraband response occurs just because 
of the change in the the population difference. Using 
(I53|, the intraband contribution of the current is 


where gs and are spin and valley degeneracy factors 
respectively. This equation can be simplified to a close 
form expression for the linear intraband conductivity 




'intray^P 


(Wp) 


QsQv ksT 

h 47r h{iujp -i- r) 




( 62 ) 


Sintra = -eVp ^ 

k 


B. Interband linear response 


then, the n’th order nonlinearity due to pure intraband 
process can be obtained using Taylor expansion: 

(59) 

k 

An intuitive symmetry argument shows that in graphene 
as a centrosymmetric crystal, even orders of nonlinearity 
do not exist. In the succeeding sections the derivation 
of linear and third order conductivity of graphene is dis¬ 
cussed. 


Interband linear optical response of graphene can be 
obtained using the master equation ((55l) . Linear optical 
response is a single photon process and unlike higher or¬ 
der terms it can be obtained independetly for interband 
and intraband contributions. The first order solution of 
the Eq. (1551) can be derived by replacing Wk and rck with 
and 0 respectively 


U'mier (k, CJp) — 


VF _ (72 -I- iujp) ^ , 

fi + {5kY Wp - 2f72a;p - 


(63) 


VI. LINEAR OPTICAL RESPONSE OF 
GRAPHENE 


Integration over the reciprocal space and including the 
density of states gives 


A. Intraband linear response 

Eq. (15^ for n = 1 gives the intraband conductivity 
tensor in the /c-space 




.(wp) — 


gsgv 


^+00 


h 47r 


df- 


(72 -I- iujp) 


u)'^ - 2ij2UJp - 


\J£‘^ (Ago — 

[f{£)-f{-£)\ (64) 


=(i) 

^intra 


(k,Wp) 


ihhjp dk 


( 60 ) 


where we have assumed that the electric field is E = 
Epexp(ia;pt) and Wp) is defined via 


^mira(k;a)p) ■ E 


It is quite obvious that in the linear regime the effect 
of spin orbit coupling can be neglected as the density 
of states on the energy axis linearly goes to zero and it 
removes the singularity and thus -I- (A^o — \rY 

can be safely replaced by unity. 

Eqs. (1621) and (IMl) are identical to the results obtained 
from linear response theory and Kubo formulation [^. 


Again the effect of radiation loss and collision in trans¬ 
port equation can be crudely incorporated in the equa¬ 
tions in a phenomenological level. We assign an imagi¬ 
nary part to the frequency, i.e. the frequency Wp will be 
replaced by ujp — iT. This imaginary portion can be deter¬ 
mined by fitting the analytical results to experimentally 
measured data. The main disadvantage of this method is, 
of course, certain lack of rigour, the method is, however 
less laborious. Performing integration in the reciprocal 
space, the off-diagonal terms vanish, we arrive at 


hntraV^P 


(Wp) = 


QsQv ^ 


1 


47r {ioJp + 

p + CO 

/ d££ 


r) 

df{£) 

d£ 


dfi-sy 

d£ 


(61) 


VII. THIRD ORDER FREQUENCY MIXING IN 
GRAPHENE 

Graphene as a centrosymmetric crystal does not ex¬ 
hibit second order nonlinearity and therefore the first 
nonlinear term is the third order one. In fact, to in¬ 
duce an optically biased second order response, transla¬ 
tion symmetry must be broken. From what has been dis¬ 
cussed so far, we know that the (/?k dependence of the n’th 
order optical conductivity tensor in the reciprocal space, 
i.e. CTk"\ appears as T„((/3k) = ctia 2 ---&n+i where di 
can be either k or (^k- It is straightforward to show that 
r„((y3k)d</3k vanishes when n is an even integer. 

Throughout this section we assume that three com¬ 
plex fields with the time dependence of and 
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g^j-g j] 2 ixing through the third order conductivity of 
graphene. As mentioned earlier, the third order opti¬ 
cal response can be interpreted as a three-photon pro¬ 
cess and three different terms may contribute in the 
third order conductivity tensor: pure intraband term 

=(3) =(3) 

f^mtra(wr,W 5 ,a;p), pure interband term w,, Wp) 

and combination of the both Wg, Wp). 


A. Intraband third order nonlinearity 

The third order intraband optical conductivity can be 
obtained using Eq. (l59ll . Like the linear response term, 
an imaginary part would be assigned to the frequency 

I 


components to include the effect of loss and collision. 
Exploiting the symmetry of ^ the leading term in 
the third order nonlinearity term is given by 

=( 3 ) 

^mira (^7 l^p) 

3\h3 {lUJp + T){i 0 Jg+r){lUJr+T) dk3 ^ > 

The notation (..., Wi-lwglwp) explicitly shows that the 
photons are ordered. This expression should be recog¬ 
nized from the third order intraband optical conductivity 
=( 3 ) 

Wr, Wg, Wp) which includes all possible permuta¬ 
tions of the incoming photons. The formal integration of 
(I65|) in the reciprocal space gives the third order intrband 
conductivity tensor 


=(3) /■ ^ _ 

^intra 


£1 J_ 


Qs9v 

487r 


(sTd + To) Vi 


[keTf 


exp 


( _ 

V fcsT 


h^iiujp + r)(za;g -|- r)(iWr + E) 


1 -I- exp ( - 


( ^bt) 


( 66 ) 


Here we have made use of the intrinsic permutation op¬ 
erator Vf whose meaning is that the right of it is to 
be summed over all possible permutations of the input 
frequencies uip, ujq and Ur- The tensors Td and To are 
defined as 

Td =xxxx-byyyy (67) 


xxyy -1- yyxx -h xyyx -h yxxy -h xyxy -b yxyx (68) 


The parameter Ec = (KbT)^/{ ehvp) is the characteristic 
electric field strength. In Eq. (l66)l . all terms inside the 
bracket are dimensionless and therefore Eo gives a rough 
estimate of the sufficient field stength to cause nonlinear 
effects. A glance shows that this coefficient is analogous 
to the Schwinger limit in QED giving the scale of the 
electric field that Maxwell’s equations are expected to 
become nonlinear. 

I 


B. Interband third order nonlinearity 

Pure interband third order nonlinearity can be ob¬ 
tained using the master equations (15511 and (|49l) in the 
moving frame. Power expansion of the inversion and the 
polarization in terms of the exciting field as Eqs. (l56)l and 

j—I 

(I57p gives Wjj.'. Assuming that the first photon Wp pro¬ 
vides time variation for (wp), the first nonzero oscilla¬ 
tory component of the Wk as well as the third harmonics 
of can be found from the following equations 

zi-k ^ + 7iw^k ^ = 4(wg)z;^^^(wp) (69) 


..( 3 ) 


,•,(3) 


2 ( 3 ) 

= 


-f 272 z}- 

- 72Ck(Wr)w|f^ - ■?k(Wr)wjf^ - 5k(Wr)Wk^^ (70) 


These equations can be solved simultaneously to find the 
third order interband conductivity tensor in the recipro¬ 
cal space which is given in Eq. (EB). Performing integra¬ 
tion in the reciprocal space and applying the permuta¬ 
tion operator yields the final expression for the interband 
conductivity tensor. 


=(3) 


.(k, Wr-( 


.)=- 


vpe 




[72 + i{uJp +UJq+ UJr)] (72 + 


[(Wp +iOq+ ujrY - 2z72(a;p + uiq+ uJr) - [z(a;p -b Wg) -b 7l] [Wp - ‘ 2 i'j 2 ^p - 


■<^k<^k<^k<^k (71) 
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=(3) 


' inter 
/*+oo 


{wr,u;„u;p)- ^ { 


h [72 + i(wp + w, + Wr)] (72 + *Wp) 


ksT [i{ujp + ujg) + 71 ] 


d£ 


[keTf 


/o + (^so - ^rY [(Wp + Wq + Wr)^ “ ‘^il 2 {^p +^q+ OJr) “ 


[/(£)-/(-£)] (72) 


The significant role of band gap opening and the band 
renormalization due to the spin-orbit coupling can be ob¬ 
served in Eq. (1721) . It is quite interesting that the integral 
given in Eq. (|72l) possesses a first order singularity in the 
absence of spin-orbit coupling. 

C. Intraband-Interband third order nonlinearity 

According to Eq. (1571) . the two level problem in the 
moving frame can provide the odd order of nonlinear 
response in terms of the effective dipole in the moving 
frame However, the motion of the frame can also be 
incorporated into the frequency mixing process. Through 
a second order intraband process the first and the second 

I 


photons whose frequencies are Wp and Wq are contributing 
to the frequency mixing which accounts for the motion of 
the primed frame. Let us assume that, the third photon 
iOr is absorbed through a first order interband process 
engendering + ujq + uir) ■ An intuitive argument 

shows that the first two photons mainly contribute in 
changing the population difference. The Taylor expan¬ 
sion of the iCk around the equilibrium state and plugging 
it into Eq. (1551) yields the third order conductivity tensor. 
The resultant equation is given in Eq. (l73l) . In Eq. (l73l) 
the higher order loss terms such as 72 r and T^ have been 
neglected. Performing integration in the reciprocal space 
and applying the permutation operator yield the final 
expression for the third order conductivity tensor. 


^intra—inter \^Q l^p) 

e^vp 72 + 2r -K i{ujp + Wg+ Ur) _yk<^kkk_ 

2tfi (iwp-I-r)(ia;g-I-r)fc dk"^ {ujp-\-WquJrY — 2i'yi{ijjp + ujq + ujY) — 


=( 3 ) 

^ intra—inter 


(cOr, 


^q; ^p) 



y 9sgv (72 + 2r) -I- i{ujp -I- a;g -I- a;,-) j, 
32tt H{ujp — iT){ojq — iV) ^ 



dS-^ [f{£) - fi-£)] 


{ujp UJq UJrY 


(ksTY _ 

2z7i (^p ^q ^r) ^k 


(74) 


where the tensor Tc is defined as 

Tc = 3(xxyy -|- yyxx) -I- xxxx 

+ yyyy - xyxy - yxyx - yxxy - xyyx (75) 


VIII. RESULTS 

The nonlinear optical behavior of graphene has been 
experimentally investigated by several groups [l^ . [^[^ . 
Four-wave mixing experiment [ 3 ^ and Z-scan technique 
have been among the most common methods to mea¬ 
sure the nonlinear response of graphene. The results of 
the experiments confirm that graphene has an exception¬ 
ally high third-order susceptibility over a wide range of 


frequency. Depending on the measurement method and 
the sample quality, various research groups have reported 
different values for the bulk susceptibility and nonlinear 
refractive index of graphene. Despite the discrepancies 
between the different measurement results, all reliable 
experiments unanimously demonstrate that the nonlin¬ 
ear response of graphene is several orders of magnitude 
stronger than that of all known semiconductors [ 1 ^ . 

The theoretical predictions for the linear and nonlinear 
optical response of graphene are shown in Figsl2] and [3l 
In our calculations, the unknown parameters are selected 
according to the experimental results. The phenomeno¬ 
logical intraband scattering rate T and two-level inter¬ 
band relaxation coefficients (71 and 72 ) obviously depend 
on the frequency, temperature and quality of the sample. 
A full theoretical investigation of the Drude conductivity 
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of graphene and possible origins of the relaxation coef¬ 
ficients is given in 13911. A ccording to the experimental 
results reported in 1^ , the relaxation coefficients are 
typically around tens of meV (r, 7 i _2 ~ lOmeV). To 
highlight the resonant features of the linear and nonlin¬ 
ear response, it is assumed that graphene is held at the 
temperature of T Ri OK. It is also assumed that the 
Fermi energy level is around lOOmeV. The energy of the 
incoming photon(s) are normalized to the Fermi energy. 

The linear optical conductivity of graphene is shown 
in Fig. [21 It is observed that the optical absorption of 
graphene is universal and independent of the frequency 
for the photon energies of fuLi ‘^Ef. The univer¬ 
sal behavior of the optical conductivity in graphene can 
be explained by the two-dimensionality and the invari¬ 
ance of the condensed matter system [l3, Sll- However, 
the dispersionless character of the absorption spectral is 
not topologically protected and even the inclusion of the 
higher order terms -such as triangular warping Q- can 
deviate it from the universally flat response. 



Tiuj/Ef 


FIG. 2. Linear optical conductivity of graphene for normally 
incident plane wave. The parameter cro = e^/(4/i) is the uni¬ 
versal optical conductivity. 


around hui ^ 2i?/. The nature of this resonant behavior 
can be explained by looking at the imaginary of the in¬ 
tegrand in Eq. (EH) . The integrand for ujp^q^r = ^, 0 ;, —w 
has two second order poles at Hk ~ ±w . The ampli¬ 
tude of the absorption for such poles is roughly propor¬ 
tional to the slop of population difference. Those reso¬ 
nances are significantly stronger around the Fermi energy 
level where the population difference abruptly changes. 
A physical interpretation of this behavior can be pre¬ 
sented based on two photon absorption of a bunch of two 
level systems in the reciprocal space. 

The equivalent third order bulk susceptibility of 
graphene is related to the third order surface dynamic 
conductivity via 


X^^\ujp,UJq,UJr) = - 


(3) ( \ 

^xxxx 1 ) 

i{uJp + UJq + UJr)dgrSo 


(76) 


where dgr is the equivalent thickness of graphene which 
is typically around d Ri SA mill and eo is the free 
space permittivity. Obviously for the case of graphene, 
the definition of the nonlinear bulk susceptibility is am¬ 
biguous due to the arbitrariness in the definition of the 
thickness of the two-dimensional structure. In the Kerr- 
type nonlinear response, the dependence of the complex 
refractive index n on the intensity of light / is given by 


n = no + n2l 


(77) 


Where / = 2eoRe{no}c |E|^ (c is the speed of light). The 
nonlinear coefficient n 2 is related to the bulk susceptibil¬ 
ity as [12] 


n 2 


4eoC I no I' 




.Im{no} 

1 — 2 - - - 7 

Re{no} 


(78) 


Fig. [3] shows the frequency dependence of the third or¬ 
der conductivity of graphene. It can be shown that the 
third order conductivity tensor is mainly dominated by 
the interband transitions and Rabi-type oscillations de¬ 
termine their frequency dependence for high-energy pho- 


i'O'N 

tons. The frequency variation of the Gxxxxioj, uj,uj) which 
is responsible for the third harmonic generation is shown 
in Fig. 3(a) The resonant features of this response can 
be explained based on the Rabi oscillations in the recip¬ 
rocal space. The integrand in Eq. EH for UJp,q,r = W 
possesses four simple poles at Hk ~ ±a;,±u;/3. In the 
absence of Pauli blocking , i.e. for fiHk 2i?y, the 
interband transitions take place. The overall response 
is the superposition of the broadened resonances in the 
reciprocal space. The resonances around the hu) ~ 2Ef 
are stronger leading to appearing of the peaks around 
fvjj ~ ^Ef and huj ^ 2,Ef. 


(o\ 

Our theoreti cal p rediction of axxxx{u!,u!,—ui) is also 
plotted in Fig. |3(b)[ This part of the nonlinearity con¬ 
tributes in the nonlinear refractive index. This compo¬ 
nent of the nonlinear response exhibits resonant behavior 


It is easy to show that for the case of graphene this ex¬ 
pression is merely independent of the particular choice 
of dgr and it introduces an intrinsic parameter. Fig. 0] 
displays our theoretical prediction for the real and imag¬ 
inary parts of the Kerr coefficient 712 at the room tem¬ 
perature. 

As the figure shows in a wide range of frequency and far 
from the resonances, the nonlinear coefficient is around 
712 ~ This results is in a reasonably good 

agreement with the experimental results provided in [38| . 
This curve also indicates that the Kerr nonlinearity of 
graphene is much stronger than all known nonlinear semi¬ 
conductors such as GaAs, Ge and AlGa. According to 
the experimental results presented in Ref. the nonlin¬ 
ear Kerr coefficient for those materials is in the order of 
which is obviously much smaller than that 
of graphene. It is quite interesting that nonlinear Kerr 
coefficient can be tuned by changing the Fermi energy 
level and extremely high Kerr type nonlinearity can be 
achieved in appropriately gated graphene monolayer and 
chiral multilayer graphene. 
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(a) 

FIG. 3. The third order conductivity of graphene normazlied 

( 3 ) ( \ 

-CO). 



flUl/Ef 


FIG. 4. The Kerr nonlinear coefficient of graphene (n 2 defined 
in Eq. GZI)) at T — 300K. Over a wide range of frequency 
the Kerr coefficient is around 10“^^ni^W“^. 


IX. CONCLUSION 

A semiclassical theory of light-graphene interaction in 
linear and nonlinear regimes has been detailed. Focus¬ 
ing on the scale-invariancy and chiral character of Bloch 
quasi-particles in the graphene Hamiltonian, the Semi¬ 
conductor Bloch Equations have been formulated. The 
advantage of SBEs is two-fold: first they provide a conve¬ 
nient mathematical scheme leading to the analytical ex¬ 
pressions for different contributions of the linear and the 
nonlinear optical response of graphene. Moreover, SBEs 
encode the topological properties of the band structure 
in an effective dipole expression appearing in the equa¬ 
tions. The mathematical structure of this effective dipole 
reveals the distinct optical response of graphene. 

The exotic transport and optical properties of 
graphene can be attributed to the chirality of the charged 



(b) 

to (To = e^/(4ft) in the unit of ^ (a) a^Jxx{ix!,u!,uj) and (b) 


carriers. Giving the mathematical description of the chi¬ 
rality and employing SBEs, it has been shown that this 
chirality leads to a remarkably strong nonlinear optical 
response. 

Using SBEs, the problem of the interaction can then 
be decomposed into the quasicalssical Boltzman trans¬ 
port and the interband time evolution. The nonlinear 
parts of the optical response can be classified as pure in¬ 
traband, pure interband and combination of the both. 
Introducing a novel mathematical framework, analytical 
expressions for different contributions of the conductiv¬ 
ity tensors have been derived for the first time. We have 
shown that a suitable band renormalization is required to 
remove a singularity in the third order interband optical 
response of graphene resulting form the chiral symmetry 
and the scale invarince of the band structure. 

The third order susceptibility and nonlinear refractive 
index have been calculated. It is shown that our predic¬ 
tion for the Kerr nonlinear coefficient for graphene is in 
reasonably good agreement with the experimentally ob¬ 
tained results. It has been demonestrated that Kerr non¬ 
linear coefficient of graphene can be tuned by changeing 
the Eermi energy level and significantly strong Kerr non¬ 
linearity can be attained in a gated graphene monolayer. 
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